###
# Replication file for everything but STM models in
# "What Makes Foreign Policy Teams Tick.."
# 5/19/2018
###

#Dependencies to Run This.
#1) data.table, ggplot2,mgcv (+ quanteda if you want to rerun linguistic data)
#2) a "scripts" folder where this lives, a "data" folder
#3) inside "data" it assumes it will find:
### fig2_public.csv
### fcast_public.csv
# the code is written assuming it is sourced to the scripts folder.

#It will generate a figures/ folder with Figures 1-8 and Figure 13

library(data.table)
library(ggplot2) 
library(mgcv)

####
#general use functions
####

#a function we use for color palettes
g_color_hue <- function(n) {
  hues = seq(15, 375, length=n+1)
  hcl(h=hues, l=65, c=100)[1:n]
}


###########
# Figure 2: Regression Discontinuity
###########
# this uses separate data so we run it first

#load- the prepared data
scores.yr1.yr2 <- fread("../data/fig2_public.csv")
#note that scores are negated so that larger values mean more accurate
ggplot(data = scores.yr1.yr2, aes(x = -centered_decision_score_yr1, y = -mean_std_score_yr2, color=treatment, shape = super_yr2f)) +
  stat_smooth(inherit.aes = FALSE, aes(x = -centered_decision_score_yr1, y = -mean_std_score_yr2, color = treatment), 
              method = "loess", span = 1, size = 1) +
  geom_point() +
  scale_shape_manual(values = c(16,1), guide=guide_legend(title="")) +
  guides(color = guide_legend(title="Treatment")) +
  theme_bw() +
  labs(x = "Year 1 Accuracy (Centered Decision Score)", y = "Year 2 Accuracy (Mean Std. Score)") + 
  coord_cartesian(xlim = c(-0.2, .085), ylim=c(-.75,.8))

ggsave("../figures/regr_discon_hypoth_1_zoomed_2_lines_reverse_score_loess.pdf", width = 7, height = 4)
#clean up the file we used.
rm(scores.yr1.yr2)


###Prep 
dat <- fread("../data/fcast_public.csv")

# Some people shifted conditions and so we drop them.
# Results are robust to leaving them in the model.
# (the nchar code just picks up the fact that a message
#  is present which indicates an exclusion of some kind)
dat <- dat[!(nchar(dat$exclude)>0 & dat$year==3),]

#flip the score so it is more interpretable
# (i.e. positive numbers are better)
dat$std_score<- -dat$std_score

#give this a better name
team_data <-dat
rm(dat)


###########
# Figure 1: Average (standardized) Brier score by group type
###########

SuperTeam<- t.test(team_data$std_score[team_data$Condition=="Top Team"])
Team <- t.test(team_data$std_score[team_data$Condition=="Team Train"|team_data$Condition=="Team No Train"])
Individual <- t.test(team_data$std_score[team_data$Condition=="Individual No Train"|team_data$Condition=="Individual Train"])

pdf(file="../figures/AccuracyByCondition-reduced.pdf", height=3.5,width=6)
par(mar=c(5,4,1,2)+0.1)
plot(Individual$estimate, 1, pch=16, xlim=c(0,.6), ylim=c(0,5), yaxt="n", ylab="Condition",
     xlab="Average (Standardized) Brier Score")
segments(Individual$conf.int[1], 1, Individual$conf.int[2], 1)
text(.15,1, "Individuals")
points(Team$estimate, 3, pch=16)
segments(Team$conf.int[1], 3, Team$conf.int[2], 3)
text(.17,3, "Team")
points(SuperTeam$estimate, 5, pch=16)
segments(SuperTeam$conf.int[1], 5, SuperTeam$conf.int[2], 5)
text(.45,5, "Top Team")
dev.off()




###########
# Figure 3: Extensive engagement: number of responses by IFP
###########

#Create a Dataset We Will Use for Several Figures
#Average by individual.
## nresponses=  Number of responses per IFP which they responded to.
## median_word_count= Median over all their replies
## median_first_word_count= Median over all their first responses
## frac_non_first_response= Fraction of responses that aren't the first
small <- team_data[,list(nresponses=nrow(.SD)/length(unique(ifp_id)), 
                         median_word_count=median(as.numeric(word_count)),
                         median_first_word_count=median(as.numeric(word_count[response.order==1])),
                         frac_non_first_response=sum(response.order>1)/nrow(.SD),
                         frac_words_after_first=median(1- sum(word_count[response.order==1])/sum(word_count)),
                         normcheck=length(unique(ifp_id))),
                   by=list(user_id,Condition)]
#Redo some of the factors to provide different condition groupings
#first a binary
small[,Condition2:=ifelse(Condition %in% c("Individual No Train", "Individual Train"), "Individual", "Team")]  
#then the more refined version
small[Condition%in%c("Individual No Train", "Individual Train"), Condition:="Individual"]
small$Condition <- factor(small$Condition, levels=c("Individual", "Team No Train", "Team Train", "Top Team"))


pdf(file="../figures/nresponsesByCondition.pdf", width=10, height=6)
ggplot(small, aes(nresponses, fill = Condition)) +
  scale_x_log10() +
  geom_density(alpha = 0.2) + xlab("Average Number of Responses per Active IFP") +
  theme_bw(base_size=18)
dev.off()


###########
# Figure 4: Median number of words used by individuals in their first response to an IFP
###########
#a few zeroes are dropped out here.
pdf(file="../figures/medianWordCountByCondition.pdf", width=10, height=6)
ggplot(small, aes(median_first_word_count, fill = Condition)) +
  scale_x_log10() + 
  geom_density(alpha = 0.2) + xlab("Median Word Count of First Response")  + 
  theme_bw(base_size=18)
dev.off()

###########
# Figure 5: Average number of explanations per active IFP by word count
###########
pdf(file="../figures/responsesWordsByCondition.pdf", width=10, height=10)
ggplot(small) + 
  geom_point(aes(y=median_first_word_count,x=nresponses)) +
  geom_density2d(aes(y=median_first_word_count,x=nresponses,colour=Condition)) + 
  facet_wrap(~Condition) + 
  scale_x_log10() + 
  scale_y_log10() + 
  xlab("Average Number of Responses per Active IFP") + 
  ylab("Median Word Count of First Response") +
  theme_bw(base_size=18) +
  theme(legend.position = "none")
dev.off()  


###########
# Figure 6: Fraction of total words written after first response.
##########

#Plot the proportion of words after the first response
#Figure 6
pdf(file="../figures/fractionOfWordsAfterFirst.pdf", width=20, height=6)
ggplot(small, aes(frac_words_after_first, fill = Condition)) +
  geom_density(alpha = 0.2) + xlab("Median Proportion of Words After First Response")  + 
  theme_bw(base_size=18) + facet_grid(~Condition2)
dev.off()

###########
# Figure 7: Ratio of responses by most prolific team member
###########
#ratio of most prolific to least
small <- team_data[Condition %in% c("Team No Train", "Team Train", "Top Team"),
             list(ratio= max(tabulate(as.factor(user_id), nbins=unique(team_size)))/mean(tabulate(as.factor(user_id), nbins=unique(team_size)))),
             by=list(ifp_id, ctt, year, Condition)]
small$Condition <- factor(small$Condition, levels = c("Top Team","Team No Train", "Team Train"))
pdf(file="../figures/ratioRespond.pdf", width=10, height=6)
ggplot(small, aes(ratio, fill = Condition)) + scale_fill_manual(values=g_color_hue(4)[c(2:4)]) + 
  geom_density(alpha = 0.2) + xlab("Ratio of Most Prolific to Average Responder")  + 
  theme_bw(base_size=18)
dev.off()

###########
# Figure 8: Extensive engagement and performance
###########

#Setting Up the Team-IFP level variables
small <- team_data[Condition %in% c("Team No Train", "Team Train", "Top Team"),
                   list(Condition=unique(Condition),fracreply=length(unique(user_id))/unique(team_size),
                        score=mean(std_score)), 
                   by=list(ifp_id, year, ctt)]
#Now aggregate to just teams
teamlevel <- small[, list(Condition=unique(Condition),
                          fracreply.median=median(fracreply),
                          fracreply.mean=mean(fracreply),
                          score=mean(score)), by=list(year, ctt)]

mod <- gam(score ~ s(fracreply.median, bs="cr") + Condition, data=teamlevel)


pdf(file="../figures/teamModelWithCondition.pdf", width=6, height=4)
par(mar=c(5,4,1,2)+0.1)
plot(mod,xlab="Median Fraction of Team Replying", ylab="Average Team Score over IFPs")
dev.off()


###########
# Appendix: Average ratio of syllables to words in explanation
###########

small <- team_data[year==2 & response.order==1 & !is.na(nsyl),]
small[,ratio:=nsyl/word_count]

Individual <- t.test(small$ratio[small$Condition%in%c("Individual Train", "Individual No Train")])
UntrainedTeam <- t.test(small$ratio[small$Condition=="Team No Train"])
TrainTeam <- t.test(small$ratio[small$Condition=="Team Train"])
SuperTeam<- t.test(small$ratio[small$Condition=="Top Team"])

pdf(file="../figures/SophisticationByCondition.pdf")
plot(Individual$estimate, 1, pch=16, xlim=c(1,2), ylim=c(0,5), yaxt="n", ylab="Condition",
     xlab="Average Syllables / Words")
segments(Individual$conf.int[1], 1, Individual$conf.int[2], 1)
text(1.25,1, "Individual")
points(UntrainedTeam$estimate, 2, pch=16)
segments(UntrainedTeam$conf.int[1], 2, UntrainedTeam$conf.int[2], 2)
text(1.25,2, "Untrained Team")
points(TrainTeam$estimate, 3, pch=16)
segments(TrainTeam$conf.int[1], 3, TrainTeam$conf.int[2], 3)
text(1.25,3, "Trained Team")
points(SuperTeam$estimate, 4, pch=16)
segments(SuperTeam$conf.int[1], 4, SuperTeam$conf.int[2], 4)
text(1.25,4, "Super Team")
dev.off()



########
# Linguistic quantities like number of syllables used above are precomputed.  
# For reproducibility we show how to recalculate them using quanteda below.
# However the code from here down can be a bit slow and isn't strictly necessary.
########
dat <- fread("../data/fcast_public.csv")
library(quanteda)
corpus <- corpus(dat$body)
#count words (remove punctuation to not count periods and commas)
#(remove symbols so errant html tags only count once and not several times
# see for example the two lines below to see how this goes wrong)
#tokens(corpus[156324])
#tokens(corpus[156324],remove_punct=TRUE, remove_symbols=TRUE)
tdata <- tokens(corpus, remove_punct=TRUE, remove_symbols=TRUE)
length(tdata)==nrow(dat) #check we haven't lost any
dat$word_count <- unlist(lapply(tdata, length))

#syllable count
nsyl <- nsyllable(tdata)
length(nsyl)==nrow(dat) #check we haven't lost any
#we want to sum these up in a way that all NAs still returns NA
dat$nsyl <- lapply(nsyl, function(x)
  if(sum(is.na(x))==length(x)) return(NA) 
  else sum(x, na.rm=TRUE))

#Flesch.Kinkaid
FK <- textstat_readability(corpus, measure="Flesch.Kincaid")
#FK does this annoying thing where it drops documents

#create an index to order the dataset
dat$order <- 1:nrow(dat)
#make a key that matches FK
dat$document <- sprintf("text%i", 1:nrow(dat))
setkey(dat, document)
#create the data table and set key for FK
FK <- data.table(FK)
colnames(FK)[2] <- "FKgrade"
setkey(FK, document)
#merge and then reorder them
dat <- FK[J(dat)]
setkey(dat, order)
